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Abstract 

Axions differ from the other cold dark matter candidates in that they form a degenerate Bose 
gas. It is shown that their huge quantum degeneracy and large correlation length cause cold dark 
matter axions to thermalize through gravitational self-interactions when the photon temperature 
reaches approximately 500 eV. When they thermalize, the axions form a Bose-Einstein condensate. 
Their thermalization occurs in a regime, herein called the 'condensed regime', where the Boltzmann 
equation is not valid because the energy dispersion of the particles is smaller than their interaction 
rate. We derive analytical expressions for the thermalization rate of particles in the condensed 
regime, and check the validity of these expressions by numerical simulation of a toy model. We 
revisit axion cosmology in light of axion Bose-Einstein condensation. It is shown that axions are 
indistinguishable from ordinary cold dark matter on all scales of observational interest, except when 
they thermalize or rethermalize. The rethermalization of axions that are about to fall in a galactic 
potential well causes them to acquire net overall rotation as they go to the lowest energy state 
consistent with the total angular momentum they acquired by tidal torquing. This phenomenon 
explains the occurrence of caustic rings of dark matter in galactic halos. We find that photons 
may reach thermal contact with axions and investigate the implications of this possibility for the 
measurements of cosmological parameters. 

PACS numbers: 95.35.+d 
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I. INTRODUCTION 



One of the outstanding problems in science today is the identity of the dark matter of 
the universe [l| . The existence of dark matter is imphed by a large number of observations, 
including the dynamics of galaxy clusters, the rotation curves of individual galaxies, the 
abundances of light elements, gravitational lensing, and the anisotropics of the cosmic mi- 
crowave background radiation. The energy density fraction of the universe in dark matter 
is observed to be 23%. The dark matter must be non-baryonic, cold and collisionless. Non- 
baryonic means that the dark matter is not made of ordinary atoms and molecules. Cold 
means that the primordial velocity dispersion of the dark matter particles is sufficiently 
small, less than about 10~^ c today, so that it may be set equal to zero as far as the forma- 
tion of large scale structure and galactic halos is concerned. Collisionless means that the 
dark matter particles have, in first approximation, only gravitational interactions. Particles 
with the required properties are referred to as 'cold dark matter' (CDM). The leading CDM 
candidates are weakly interacting massive particles (WIMPs) with mass in the 100 GeV 
range, axions with mass in the 10~^ eV range, and sterile neutrinos with mass in the keV 
range. To try and tell these candidates apart on the basis of observation is a tantalizing 
quest. 

For a long time, it was thought that axions and the other forms of cold dark matter 
behave in the same way on astronomical scales and are therefore indistinguishable by obser- 
vation, whether it be observations of large scale structure or measurements of cosmological 
parameters. More recently, however, it was pointed out that dark matter axions form a 
Bose- Einstein condensate (BEG), as a result of their gravitational self- interactions, when 
the photon temperature reaches about 500 eV 0]. Axions or axion-like particles are special 
because they are a degenerate Bose gas. The other dark matter candidates, which we refer 
to henceforth as 'ordinary cold dark matter', are non-degenerate. This raises the question 
whether axions may be observably different after all. It was shown in ref. 0] that, on 
all scales of observational interest, density perturbations in axion BEG behave in exactly 
the same way as those in ordinary cold dark matter provided the density perturbations are 
within the horizon and in the linear regime. On the other hand, when density perturbations 
enter the horizon, or in second order of perturbation theory, axions generally behave differ- 
ently from ordinary cold dark matter because the axions rethermalize to let the axion state 
(i.e. the state most axions are in) track the lowest energy state 0]. 

Axion Bose-Einstein condensation appears to resolve a puzzle that has arisen in the study 
of the inner caustics of galactic halos. The structure of the inner caustics depends on the 
angular momentum distribution of the infalling particles. If the particles fall in with net 
overall rotation, the inner caustics are rings whose cross-section is a section of the elliptic 
umbilic (-D-4) catastrophe, called 'caustic rings' for short 0, 3. If the velocity field of 
the infalling particles is irrotational, the inner caustics have a 'tent-like' structure which is 
described in detail in ref. jsj and which is quite distinct from that of caustic rings. The 
radii of the caustic rings, assuming that the dark matter falls in with net overall rotation, 
were predicted [3] using the self-similar infall model of galactic halos [6|, generalized to allow 
angular momentum for the infalling particles 0, [sj . Evidence was found for the existence 
of caustic rings at the predicted radii. The evidence is summarized in ref. jsf. Now, the 
puzzle is that ordinary cold dark matter has an irrotational velocity field [sl and is therefore 
incompatible with the existence of caustic rings. Axion Bose-Einstein condensation resolves 
the puzzle provided the axions rethermalize sufficiently quickly that most of them go to the 
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lowest energy available state before falling in [2|, |9[. The lowest energy state consistent with 
the total angular momentum the axions will have acquired from neighboring inhomogeneities 
by tidal torquing is a state of rigid rotation on the turnaround sphere, the simplest form 
of net overall rotation. Therefore, if the dark matter is an axion BEC that rethermalizes 
sufficiently quickly, the inner caustics are rings. Furthermore it was shown jol that the 
caustic rings are all in the same plane (that of the galactic disk), that the overall size of the 
rings is predicted correctly by tidal torque theory, and that the relative sizes of the rings are 
precisely as predicted by the self-similar model and therefore consistent with the evidence 
for caustic rings published earlier. One thing remains to be done: to show that axions about 
to fall into a galactic potential well rethermalize sufficiently quickly that they almost all go 
to the lowest energy available state. It is one of the goals of this paper to understand cosmic 
axion thermalization sufficiently well to be able to verify this. 

There is a second motivation for studying the thermalization rates of dark matter axions in 
detail 10|]. We mentioned that cold axions thermalize by their gravitational interactions. But 
gravity is universal. If axions thermalize by gravitational interactions, they may also enter 
into thermal contact with other species. If the axions enter into thermal contact with the 
cosmic photons, the photons will cool and some cosmological parameters, in particular the 
baryon to photon ratio at primordial nucleosynthesis and the effective number of neutrinos 
(a measure of the radiation density at the time of decoupling), will be modified compared 
to their values in the standard cosmological model. This opens the possibility of being able 
to distinguish axions from the other forms of cold dark matter by measuring cosmological 
parameters. 

Due to their unusual properties, dark matter axions occupy a unique and uncharted 
region in the physics of many-body systems. On the one hand, they are highly condensed, 
which greatly exaggerates the quantum effect of Bose-enhancement in scattering processes. 
On the other, the fact that their energy dispersion is very small implies that they are outside 
the realm of the 'particle kinetic regime'. In this case, the picture of instantaneous collisions 
breaks down, so that the usual Boltzmann equation no longer applies. Together, these 
properties - high occupation number and small energy dispersion - are highly atypical, and 
thus very little attention has been paid to the study of such systems. Existing techniques 
in non-equilibrium statistical mechanics are not applicable to dark matter axions, rendering 
the estimation of their thermalization rate non-trivial. 

The outline of this paper is as follows. In Section II, we give a definition of cold dark 
matter and ask under what conditions axions behave as such. We find that axions behave as 
CDM on all scales of observational interest except when they thermalize or rethermalize. In 
Section III, we derive expressions for the thermalization rate of particles in the 'condensed 
regime', i.e. when their energy dispersion is small compared to their interaction rate. The 
condensed regime is the one relevant to cosmic axion thermalization. In Section IV, we check 
the validity of our estimates of thermalization rates in the condensed regime by numerical 
simulation of a toy model. In Section V, we revisit axion cosmology in light of axion 
thermalization and Bose-Einstein condensation. Section VI provides a summary. 



II. AXIONS ARE DIFFERENT 

Both axions and WIMPs are considered forms of cold dark matter. Furthermore, until 
recently, axions and WIMPs were thought to be indistinguishable on observational grounds, 
i.e. indistinguishable on the basis of purely astronomical data. The discovery [2| that dark 
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matter axions form a Bose-Einstein condensate has changed this view since axion BEC is 
claimed to have observable consequences jl, [o], E^. This raises questions. First, precisely 
under what conditions are axions and WIMPs cold dark matter? Second, if both axions and 
WIMPs are CDM, how do they differ? There is a preliminary question: how is cold dark 
matter defined? 

The purpose of this section is to discuss these three questions in a general way, to set the 
stage for the more detailed calculations that follow. For the sake of brevity, we call WIMPs 
all cold dark matter candidates which do not form a degenerate Bose gas. As far as we are 
aware, this includes all cold dark matter candidates except axions and axion-like particles. 
In particular, sterile neutrinos are called WIMPs here. (Neutrinos are not WIMPs since 
they are hot dark matter.) 

To start off it is worth emphasizing that, at the fundamental level, axions and WIMPs 
are very different. The surprise is really that they have similar properties as far as large 
scale structure is concerned. Both axions and WIMPs are described by quantum fields. 
Furthermore, both are excellently described by classical limits of quantum fields. But the 
classical limits are different in the two cases: WIMPs are in the classical particle limit 
whereas (decoupled) axions are in the classical field limit. In the classical particle limit one 
takes h ^ while keeping E = huj and p = hk fixed. Since u,k — > oo, the wave nature of 
the quanta disappears. WIMPs are to excellent approximation classical point particles. In 
the classical field limit, on the other hand, one takes h ^ for constant u and k. E = Mhuj 
and p = Mhk are held fixed by letting the quantum state occupation number A/" — > oo. 
This is the limit in which quantum electrodynamics becomes classical electrodynamics. It 
is the appropriate limit for (decoupled) cold dark matter axions because they are a highly 
degenerate Bose gas. The axion states that are occupied have huge occupation numbers, 
M ~ 10^^ [2]. The need to restrict to decoupled axions will be explained shortly. 



So axions and WIMPs are fundamentally different even if it turns out that both can 
legitimately be called CDM. The distinction is not just academic, and is certainly important 
if axions thermalize, i.e. if axions find a state of larger entropy through self-interactions. 
Recall that, whereas statistical mechanics makes sense of the behaviour of large aggregates 
of classical particles (it was invented by Boltzmann to derive the properties of atoms in the 
gaseous state) it fails to make sense of classical fields. In thermal equilibrium every mode 
of a classical field would have average energy k-oT. As Rayleigh pointed out, the energy 
density is infinite then at finite temperature due to the contributions from short wavelength 
modes. Thus the application of statistical mechanics to classical field theory (classical elec- 
trodynamics in particular) is in direct disagreement with observation. As is well-known, the 
disagreement is removed because of, and only because of, quantum mechanics. 

What do we learn from this? If the axions are decoupled (i.e. do not interact and hence 
do not thermalize), they behave to excellent approximation like classical fields. As such, they 
are quite different from WIMPs since WIMPs are classical particles. So the argument why 
decoupled axions behave as CDM must be different from the argument why WIMPs behave 
as CDM. If the axions thermalize, they are not even described by classical fields. Instead 
they form a Bose-Einstein condensate, an essentially quantum-mechanical phenomenon. If 
the axions form a BEC, the argument why the axion BEC behaves as CDM must be different 
again. 

The main purpose of this paper is, in fact, to show that axions thermalize as a result of 
their gravitational interactions, that they form a BEC, and that a rethermalizing axion BEC 
behaves differently from CDM. This will be done in Sections III, IV and V. In the remainder 
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of the present section, we first adopt a definition of cold dark matter, next we discuss under 
what conditions WIMPs behave as CDM, and then discuss under what conditions axions 
behave as CDM. 



A. Cold dark matter 

Astronomical data, in particular data on the cosmic microwave background anisotropies 
and on the large scale structure of the universe in the linear regime, imply the existence of 
a new kind of stuff called 'cold dark matter' What explains the data may technically 
be called a perfect fluid with zero pressure The state of such a fluid is characterized 
by a density p{r,t) and a velocity field v{f,t). These four degrees of freedom satisfy the 
continuity equation 

dtp + V-{pv) = (2.1) 
and, in the Newtonian limit of gravity, the equation of motion 

dtV + {v-V)v = -V^ . (2.2) 

$ is the gravitational potential. All forces on the CDM fluid other than gravity are assumed 
to be negligible. The CDM fluid is a source for the gravitational potential: 

V2 $(f, t) = 47rG(p(f, t) + ...) (2.3) 

where the dots represent other sources. The stress-energy-momentum is given by 

7-«o = p , T'^ = pv^ , r'=' = pt;V . (2.4) 

We will take the above to be the defining properties of CDM because stuff with those proper- 



ties explains the data mentioned [Uj. We will not concern ourselves with the very interesting 
question to what extent the data demand CDM to be pressureless and collisionless. Instead 
we want to ask in what limit WIMPs and axions have the properties stated above. 

To show that WIMPs or axions behave as CDM it is sufficient to show that, in some 
average sense, they have the same stress-energy-momentum tensor, Eq. (12.41) . Indeed the 
conservation law D^T^^ = is equivalent in the Newtonian limit to Eqs. (12. ip and (12. 2p 
whereas Einstein's equations imply Eq. (12.30 . 

For the sake of clarity, let us emphasize that the above is only a good definition of 
CDM in the linear reg ime of structure formation. In the non-linear regime CDM produces 
discrete flows [H, ij and caustics [l, 0| and these cannot be described by Eqs. (12.11 - 
12. 3p . Instead the state of CDM in the non-linear regime is described by the embedding of 
a 3-dim. hypersurface in 6-dim. phase space plus the density of particles on this 3-dim. 
hypersurface. This more general (and hence better) description of CDM defines it as a 
cold collisionless fluid. In the linear regime, the two descriptions are equivalent because the 
3-dim. hypersurface covers physical space only once. 



B. WIMPs 

WIMP dark matter is a (huge) collection of classical point particles. Let their number 
be N . Their state is given at time t by giving all the particle positions ri{t) and velocities 
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Vi(t), i = 1,2, ... A^. Their time evolution from appropriate initial conditions is determined 
by the equations of motion: 

^ = -^ = -Vt(r.,<) . (2.6) 
The stress-energy-momentum tensor is 

N 

j.ok j.ki^^^^ ^) ^ , v^{t) , v^it^tMr- fj{t)) , (2.6) 

where niw is the WIMP mass. The indices i and j label the particles whereas k, I = 1,2,3 
label the spatial directions. We want to derive the conditions under which the WIMP parti- 
cles behave as CDM, defined in the previous subsection. One goes from the one description 
to the other by averaging over a volume V centered at r: 

N 

p{r,t) = V'^^[ dhT'\r + s,t) = ^ [ d'sT6{f+s-f,{t)) 
V{f,t) = kr'' = - dhkT'^\f+s,t)^'^ / dhT] vj{t) S{f+s-rj{t)) 

^ Jv y Jv 

N 

T'\r, t)^^j^ dh T^\r + I t)^^l^ d's v]{t) v\{t) 5(f + vM (2-7) 

We take V to be time independent. V must be large enough that the fluctuations in particle 
number inside V are negligible. Conservation of the averaged stress-energy-momentum T^^ 
follows merely from the conservation of the original stress-energy-momentum tensor T^'^ . 
But T'^'^ does not generally have the perfect fluid form. The velocity of the WIMP fluid is 

vm^^vm . (2.8) 

p{r,t) 

For each particle inside volume V , let us define 

l^v,{T,t)=Vi{t)-v{r,t) . (2.9) 

We have then 



T^\f,t) = p{r,t) v\f,t) v\r,t) + 



N 

ITLw 



V 



f d'^sS^ Av^^ l^v\5{r + s-rj{t)) . (2.10) 



If the velocity distribution inside V is isotropic, we may define the pressure p(f, t) by 

J d'sY, ^^j Av'jKr + s- fj{t)) = p{f, t)6''' (2.11) 

so that 

r'^' = p(f, t) v^{r, t) v\r, t) + p{r, t)5^^ . (2.12) 
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The stress-energy-momentum tensor is then said to have the perfect fluid form. If, in 
addition, there is a relation p{p) that determines pressure in terms of density (e.g. the 
relation p(p) for atomic gases implied by adiabatic expansion/compression), the conservation 
laws Dfj^T^'^ = 0, being four equations for the four unknowns p{r,t) and v{f,t), allow us to 
determine the evolution from arbitrary initial conditions. 

As mentioned, CDM is defined to be a pressureless perfect fluid. The pressure is a measure 
of the velocity dispersion of the particles: p = | < (Av)^ >. For WIMPs to be CDM, their 
velocity dispersion must be sufficiently small. The best limit comes from the fact that WIMP 
density perturbations on scales less than their free streaming distance are erased. Too large 
a velocity dispersion is inconsistent with the existence of the smallest observed large scale 
structure, that which gives rise to the 'Lyman-alpha forest'. This constrains the primordial 
WIMP velocity dispersion to be less than of order 10~®c today. If the WIMP particles have 
ordinary weak interactions, their kinetic energy decouples in the early universe [i.e. they 
stop colliding with other particles in the primordial soup) when the photon temperature is 
a few MeV. Their velocity dispersion today is then 



The lightest allowed mass is therefore of order keV, which is that of sterile neutrinos [l^. 

In summary, WIMPs behave as CDM if there exists a volume V large enough so that the 
fluctuations in particle number inside the volume are negligible, and small enough that the 
velocity dispersion of the particles inside the volume is negligible. Under these conditions, 
the WIMPs behave as CDM for any observer who does not resolve length scales smaller 
than Vs. 

Finally, let us remark that gravitational interactions among point particles is a source of 
velocity dispersion: 



where dw is the average interparticle distance. It is not possible for the velocity dispersion 
to be less than Av\g. Furthermore Av\g grows by gravitational instability as the particles 
aggregate into clumps of ever increasing size. Av\g is completely negligible for all proposed 
WIMP candidates, such as neutralinos and sterile neutrinos. However in present numerical 
simulations of structure formation the particle mass is 10"^ solar masses or more and hence 
Av\g is larger than 4 • 10~^, which is more than the allowed velocity dispersion of CDM. 
By this criterion and others [13], present simulations do not have sufficient resolution to 
describe CDM. For the particles in the simulations to have velocity dispersion as low as that 
required of CDM, the particle mass should be no more than approximately one solar mass. 

C. Axions 

In the first part of this subsection, we give the argument why decoupled cold axions 
behave as CDM. Cold axions remain decoupled in the early universe from shortly after 
their first appearance during the QCD phase transition till the time when gravitational self- 
interactions cause them to form a BEC. The Bose-Einstein condensation of axions occurs 




(2.13) 




(2.14) 
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when the photon temperature is approximately 500 eV ( ^012 Qev ) ^ (^^^ Section V). After 
the condensation is complete, almost all axions are in the same state. In the second part of 
this subsection, we give the argument why axions behave as CDM when they are all in the 
same state and that state does not change by rethermalization. 



The results obtained in this subsection are not new, although some o 
be. Descriptions of CDM in terms of a classical field are given in refs. 



of CDM in terms of the Schrodinger equation are given in refs. |20|, ff^ 



the derivations may 



16Hl9[. Descriptions 



1. Decoupled axions 

As mentioned earlier, decoupled axions behave as classical fields because the quantum 
state occupation numbers of those states that are occupied are huge, of order 10^^ [2|. In 
flat space-time, the classical axion field 0(r, t) satisfies the equation of motion 

(9^2 - + m^)^ = (2.15) 

where m is the axion mass. Space-time curvature does not play a role in this discusssion, so 
we ignore it. The stress-energy- momentum tensor is 

TOO = i ^ (V0)2 + mV^) 

T'^' = di(l) + ^ (0' - (V0)2 - m202j _ (2.16) 

To find out under what conditions decoupled axions behave as CDM, we average again over 
a spatial volume V located at r: 

r>'-{f,t) = ^[d'sT>'''{f+s,t) , (2.17) 

where /i, = 0, 1, 2, 3. To perform the average, we Fourier transfom (p within the volume V: 

0(x , t) = J2 t)e'P-^ + (P\p, r, t)e-^^-^) , (2.18) 

p 

where p ■ x = —p^t + p ■ x and p^ = + \/ {p}'^ + w?. The Fourier components depend on the 
position r of the volume. One finds: 

p(f,t) = r°° = X; N{p,r,t)p' 
p 

P(f, t) = k-p"" = J2 r, t)p 

p 

r"' = J2 N{p,f,t)p'p^\ , (2.19) 
p 

where 

iV(f,p,t)=2p°|0(p,f,t)|2 . (2.20) 
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To obain the expression for T^^ it is necessary to average not only over the volume V , but 
also over time. The time interval to be averaged over should be much longer than m~^. 
Eqs. (I2A9|) show that 

M{rJ,t) = ^^N{f,p,t) (2.21) 

may be thought of as the phase space density of the coarse grained axion field. 

We assume the axions to be non-relativistic. The average momentum at r is 
(p) {f, t) = ^V. The classical axion field behaves like CDM provided the momentum dis- 
tribution N{f,p,t) is narrowly peaked around t) for all r. The velocity field is then 

T''\f, t) = p{f, t)v\r, t)v\f, t) + y2 N{f,p, t)-6pk{f,p, t)6pi{f,p, t) (2.22) 

V 

where 

5p(f,p,t)=p-(p)(f,t) . (2.23) 

Therefore one condition for decoupled axions to behave as CDM is that they have sufficiently 
small velocity dispersion. That condition is satisfied for the axions produced during the QCD 
phase transition. Indeed their velocity dispersion is [see Section V.A] 

bv{t) V¥ 2.24 

^ ' mti a{t) ^ ' 

where ti ~ 2 ■ 10~^sec ( iqi^'^ cev ) ^ time when the axion mass effectively turns on 

during the QCD phase transition, and a(t) is the cosmological scale factor. If these axions 
remained decoupled afterward, their velocity dispersion today would be approximately 5 • 

5 

10^^^ c ( ioi2 Gev ) certainly small enough to be called CDM. 

The second condition for decoupled axions to behave as CDM is that they be observed 
only on length scales larger than their correlation length. The low velocity dispersion of cold 
decoupled axions implies a correspondingly large correlation length 

t{t) ^ t '^^^^ (2 25) 

m6v(t) a(ti) 

On scales less than i{t), axions behave differently from CDM. Indeed for CDM, p{r^t) and 
p{f ',t) are independent variables no matter how close r and r ', whereas for decoupled 
axions, p{r,t) and p{f',t) are independent variables only if |r — r '| > i. However, even 

today, i is only about lO^^cm ( iQia'^ oev ) ^ ' ^luch smaller than any scale on which we have 
observational information on the nature of CDM. 



2. All axions in the same quantum state 



When all the axions are and remain in the same quantum state, the axion fluid also 
behaves as CDM on all scales of observational interest. The argument for this was given in 
ref. j2|, but we restate it here for the sake of completeness. Note there is no requirement 
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that the axion state is the lowest energy state. In principle the axion state may be any state. 
It is important, however, that the axions remain in the same state all the time. 
The quantum axion field may be expanded in modes labeled a: 



J2 [«a + a)j (2.26) 



where the ^d{x) are the positive frequency c- number solutions of the Heisenberg equation 
of motion for the axion field 

D^^D.^ix) = g^'ld^d, - r^^^dxMx) = m\{x) , (2.27) 

and the aa and at are annihilation and creation operators satisfying canonical commutation 
relations. Eqs. (12.271) are written in curved space-time to emphasize that the modes depend 
on the background. The words 'mode' and 'particle state' are equivalent. Let us assume 
that, except for a tiny fraction, all axions go to a single particle state with label c5 = 0. The 
corresponding ^o{x) is the axion wavefunction. For the sake of brevity, we call the state of 
the axion fiuid with almost all axions in a single particle state an axion BEC whether or not 
that particle state is the lowest energy state. 

Ignoring the small fraction of axions that are not condensed into mode a = 0, the state 
of the axion field is \N >= (l/TM) (aj)^|0 > where |0 > is the empty state, defined 
by |0 > = for all a, and is the number of axions. The expectation value of the 
stress-energy-momentum tensor is 

%u{f,t) =< N\T^,\N >= N[d,%d,% + d,%d^%+g,,{-dxnd^%-m^<^*M] . (2.28) 

To see under what conditions the axion BEC behaves as CDM, we may restrict ourselves to 
Minkowski space-time. Since the axions are non-relativistic, ^o{x) = e~*'^*\E'(x) with \E'(x) 
slowly varying. Neglecting terms of order ^ dt compared to terms of order one, Eq. (I2.27p 
becomes the Schrodinger equation: 

idt^ = ^ . (2.29) 

2m 



The wavefunction may be written as [15 



^(x,t) = ^^=B(x,t)e'^^'^''^ . (2.30) 

In terms of B{x, t) and /3{x, t) the energy and momentum densities are Too = P = rn {B{x,t))'^ 
and Toj = —pvj = — {B{x,t))'^ djP, in the non-relativistic limit. The velocity field is there- 
fore v{x,t) = ^V/3(x, t) Eq. f l2.29p implies the continuity equation, Eq. (12. ip . and the 
equation of motion 

dtv'' + v^djv^ = -Vq (2.31) 

where 

?(^,i) = -^^ ■ (2.32) 
2m^ yT? 

Following the motion, the stress tensor is 

Tjk = pVjVk + '^^^^jP^'^P ~ ^jk^'^P) ■ (2-33) 
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Comparison with Eqs. (12. 2p and (12. 4p shows that axion BEC differs from CDM: the last 
terms on the RHS of Eqs. (I2.3ip and (I2.33P are absent in the CDM case. 

The extra terms are due to the Heisenberg uncertainty principle. When one attempts to 
localize the axion BEC within a region of size D, the axions acquire a minimum momentum 
spread of order h/D and hence a velocity dispersion Af ~ h/mD. The axion BEC 

tends therefore to delocalize. The tendency to delocalize is described by the force per unit 
mass — Vq appearing in Eq. (I2.3ip . and by the extra stresses in Eq. (I2.33p . CDM has no 
such tendency to delocalize. Whether this difference between axion BEC and CDM is of 
observational relevance depends on the axion mass m. The properties of CDM have been 
observed only on very large scales, D > 100 kpc. The associated velocity dispersion is tiny, 

of order 3 ■ 10^^^ c ( iL_L£Y] for the axion masses of interest to us. Hence there is no 

distinction between axion BEC and CDM, on scales of observational interest, which follows 
from the extra delocalizing forces on the axion BEC. As an illustration, consider the equation 
for the evolution of axion BEC density perturbations in the matter dominated phase of the 
expanding universe [l8| jl, [l9| : 



d^S + 2Hd,S-(^4nGpo-^y = . (2.34) 

Here, H = a/a is the Hubble expansion rate, a{t) the scale factor, po the average density of 
the axion BEC (assumed to be the dominant form of matter), k the co-moving wavevector, 
and S the Fourier component of (p — Po)/Po with wavevector k. The last term in Eq. (I2.34p 
is absent for CDM. As a result of this term, whereas the Jeans length of CDM vanishes, 
axion BEC has Jeans length [Tsl |i,[l9|: 

kj' = (16.Gp™^)-i = 1.02 . 10" cm * ('°'";^""') ' . (2.35) 

It is small compared to the smallest scales (~ 100 kpc) for which we have observations on 
the behavior of CDM. 

Several authors have proposed [2l| that the dark matter is a Bose-Einstein condensate of 
particles with mass of order 10^^^ eV. When the mass is that small, the dark matter BEC 
behaves differently from CDM on scales of observational interest as a result of the tendency 
of the BEC to delocalize. We are not considering this interesting possibility here because 
axion masses are not expected to be so small. 

In summary, cold axions with mass in the 10~^ range (give or take a few orders of 
magnitude) behave as CDM on all scales of observational interest in two different cases: 1) 
when they are decoupled and therefore behave as a classical field, and 2) when they are all in 
the same particle state. The latter result constrains the conditions under which axion BEC 
may differ from CDM. For the axion BEC to differ from CDM it must be rethermalizing, 
i.e. the state that most axions are in must be changing in time, as when it tracks the lowest 
energy available state. 



III. AXION FIELD DYNAMICS 

We consider two types of processes through which dark matter axions may thermalize in 
the early universe: self-interactions of the type and gravitational self-interactions. We 
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will see in Section V, confirming the results of ref. j2|, that the interaction is barely 
effective in thermalizing axions for a brief period just after they are produced during the QCD 
phase transition, whereas gravitational self-interactions are clearly effective in thermalizing 

axions after the photon temperature reaches approximately 500 eV ( 1012'^ cev ) ^ • When the 
axions thermalize, they form a Bose-Einstein condensate. 

In this Section, we first express the axion field dynamics in terms of a set of coupled 
oscillators. Next we obtain the evolution equations for the oscillator occupation numbers 
up to second order in perturbation theory. We show that the second order terms yield 
the usual Boltzmann equation for the evolution of the momentum distribution provided 
that the transition rate between momentum states is small compared to their spread in 
energy: F << 5uj. This latter condition defines the 'particle kinetic regime'. In the particle 
kinetic regime, the first order terms in the evolution equations are irrelevant because they 
average out in time. Dark matter axions are, however, in the opposite regime: 5u: « F, 
which we refer to as the 'condensed regime'. In the condensed regime, the first order terms 
do not average out in time and dominate over the second order terms. Using the first 
order equations, we give expressions estimating the relaxation rate of the axion dark matter 
momentum distribution due to A</)^ self-interactions and to gravitational self-interactions. 



IS 



A. Axion interactions 

Including lowest order self-interactions but neglecting gravity, the action density of axions 

Ca = -\d,cj>d^(t^-\m'(t^' + ^<P' + ... . (3.1) 

The dots represent interactions of the axion with other particles and axion self-interactions 
which are hi gher order in an expansion in powers of 0. For the axion that solves the strong 



CP problem |22|-|24l| . the mass m and self-coupling A are given by 



m^f^ yjm^ rrid „ ^^^g 10^^ GeV 
m = — ~ 6 ■ 10 eV 

fa TUu + TUd fa 

A = ^ ,""'+'"1, ^ 0,35 ^ (3.2) 

where is the pion mass, ~ 93 MeV the pion decay constant, and m„ and are the 
up and down quark masses. The formula for the axion mass [23] is obtained by expanding 
the effective potential for pions and axions to second order in the physical axion field. To 
obtain A, simply expand to fourth order. 

We introduce a cubic box of volume V = Lr" with periodic boundary conditions at its 
surface. Inside the box, the axion field (p{x,t) and its canonical conjugate field 7i{x,t) are 
expanded into Fourier components 

n 

7i{x,t) = J^i-iUri) (anit) M^) - *^(^)) (3-3) 

n 

where 

<l>ri(x) = , ^ e*P"-^ , (3.4) 
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and n = {rii, n2, n-^) with rij. (k = 1, 2, 3) integers, = ^n, and uj = a/p ■ p + rn'^- The a^j 
and at satisfy canonical equal-time commutation relations: 

[aH{t),al,{t)] = 6H,n' , [af{{t),afi'{t)] = . (3.5) 

Note that we are quantizing in the Heisenberg picture, not the interacting picture. 
Provided the axions are non-relativistic, the Hamiltonian is 



H= ^ ^ Ur^, aUri, + Yl 7 4i42«n3«n4 (3-6) 



n ni,n2,n3,n4 



where 



ni,n2 ~ 4^2y ""l+"2,n3+n4 • l-J- ' J 

The presence of the Kronecker symbol 5,i^+^2,n3+n4 expresses momentum conservation for 
each individual interaction. In Eq. f l3.6p we dropped all terms of the form a'^a'^a^a^ , a a\ 
a a a a, and a a a . We are justified in doing so because energy conservation allows only 
axion number conserving processes at tree level. Axion number violating processes occur in 
loop diagrams but can be safely ignored because they are higher order in an expansion in 
powers of 4-. In fact, all axion number violating processes, including the axion decay to two 

Jo. 

photons, occur on time scales much longer than the age of the universe in the axion mass 
range (10~^ eV) of interest. 

In the Newtonian limit, the gravitational self- interactions of the axion fluid are described 

by 

H, = -^ [ d'x pAMp($A (3.8) 
^ 2 J \x-x'\ ^ ^ 

where p = ^(vr^ + m?(fp') is the axion energy density. Because we neglect general relativis- 
tic corrections, our conclusions are applicable only for processes well within the horizon. 
Substituting and vr by their expansions in terms of creation and annihilation operators, 
Eqs. (13. 3p . and dropping again all axion number violating terms, Eq. (13. 8p becomes 

H,= E i ^ S;S 4,4.«r.3«n^4 (3.9) 

rii,n2,n3,n4 

where 

Having expressed the dynamics of the axion field in terms of coupled oscillators, we now 
study the time evolution of such systems. 



B. Evolution equations 

We just saw that the axion field is equivalent to a large number M of coupled oscillators 
with Hamiltonian of the form 



^ 1 
^ = Y1 ^A^i + 5Z 4 ^'ki 4«/«i«i ■ (3-11) 



j=l 
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In particular, the total number of quanta 

M 

N = "^a]aj (3.12) 

i=i 

is conserved. In Eq. (13. lip . A^-'^ = AH = A-j^i = A^j *. The question of interest now is the 
following: starting with an arbitrary initial state, how quickly will the averages (A4) of 
the oscillator occupation numbers A4 = a^a^ approach a thermal distribution? The usual 
approach to this question uses the Boltzmann equation. However, we will see that the 
assumptions underlying the Boltzmann equation are not valid for the cold axion fluid. So 
we need a more general approach. 

It is instructive to start with a system of just four oscillators (M = 4) and one interaction 
between them: 

4 

H = ''^^cujOjaj + A^alala^a^ + ala\aia2) ■ (3.13) 
i=i 

We have in this case 

di = i[H, ai] = i{—ujiai — Aala3a4) (3-14) 

and therefore 

A/i = iA{aia2ala\ — alala^a^j (3.15) 

and similar equations for the other aj and Afj. We solve the equations perturbatively up to 
C(A2). Let us define 

aj{t) = {Aj + B^{t))e-''^^' + ^(A^) (3.16) 

where Aj = aj{0) and Bj(t) are respectively zeroth and first order, and -Bj(O) = 0. Eq. fl3.14p 
implies 

Bi = -iAAlAsA^e^'^' + ^(A^) , (3.17) 
with Q = ui + UJ2 — u}^ — UJ4, and therefore 

9 Of 

B,{t) = -zAAlA.A, e+*™/2 - sin(^) + ^(A^) . (3.18) 

Substituting this into Eq. f l3.15p . we have 
ATi = iA{A^A2AlAle''^' - h.c.) 

+ a\aIa2A:,aIa^a\ + ^14^34^4^ 

- A^A\A2A\A^A\- A^A\A2A\A\A^)e-''''l'' - sin(y) + /i.e.] + ^(A^) .(3.19) 

Eq. fl3.19p may be recast in the form 

M = iA{AxA2A\A\e~'^^ - h.c.) 

+ A2 [AfMAfi + l)(A/'2 + 1) - MiM2{M^ + l)(A/'4 + 1)] | sin(fit) + 0(A=') (3.20) 

by rewriting the second order terms. 
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We now generalize to a system with an arbitrarily large number M of coupled oscillators, 
Eqs. (13. lip . The calculation is essentially the same as for the M = 4 toy model, except that 
there is a multiplicity of interaction terms to keep track of. One finds (/ = 1...M) 



i,j,k=l 

1 2 

+ E ^\^^^\'[M^fJiM + l){^fk + l)-M^fk{^f^ + l){^fJ + l)]^sm^^^ 

fe,i,i=i «i 

M M ^ ^ C)mn 



k,i,j = l p,m,n=l Ip 
{p;m,n)/{fc;i,j) 

k,i,j = l p,m,n = l kp 



M M 



- E E [^Aii,Ar4^l4^"^^"^.^^^"^'^^^"^'^* 7^ M^t) + h.c 



ley in ip I K p "<■ 'I- J Qjp 

(p;m,n)^(j;l,k) 



M M 



- E E [^A;iA-MMt4A.A^A.e^(^^'+^^^"/2)* -1^ sin(^t) + h.c] 

k,i,j = l P,m,n=l ^ iZmn ^ 

(p-m,n)^(i;l ,k) 

+ OiA^) , (3.21) 

where fifj = cu^ + cu; — cjj — Wj. The double sums are absent in the toy model because there is 
only one interaction in that case. At any rate, the double sums will not play an important 
role in the discussion that follows. 



C. The particle kinetic regime 

The interaction terms cause transitions in which a quantum is moved from oscillator j to 
oscillator k and at the same time a quantum is moved from oscillator i to oscillator /. As 
a general rule, energy is not conserved in any one transition. Only the total energy of the 
fluid is conserved. However, for the oscillator couplings of the axion field case, Eqs. (13. 7p 
and (I3.10p . three-momentum is conserved in each transition. 

In almost all systems of physical interest, the rate at which the occupation number of a 
typical system oscillator changes is small compared to the energy exchanged in the transitions 
it makes, i.e. fl^jt » 1. This condition defines the 'particle kinetic regime'. In the particle 
kinetic regime, the first order terms in Eq. (I3.2ip . as well as the second order terms in the 
double sums, average to zero in time. In addition, energy is conserved in each transition 
because 

-^sin(^]gt)^27^<5(^]g) (3.22) 

ij 
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for Qf^ t oo. We have then 



M 

<Afi>=+ Yl ^\^-i\^MAfj{M + l){Afk + l)-MAfk{M + l){Afj + l)]2n6{n'^j^^ 



i,j,k=l 

(3.23) 

Note that the average on the LHS of this equation is a time average, not a quantum- 
mechanical average. Eq. f l3.23p is vahd as an operator statement. 

If we substitute for the oscillator couplings those, Eq. (13. 7p . implied by self- 
interactions and take the infinite volume limit, we recover the Boltzmann equation [25| 
as an operator statement: 

y 1 f d?p2 f d^ps f d^p. 



2uJi J {2Tif2uj2 J (27r)32a;3 J (27r)32w4 
{2nY5\p,+p2-p,-p,) \ 

\{Mx + \){M2 + V)M^Ma - MxM2{Mz + l)(A4 + 1)] (3.24) 

including the factor (explicitly shown) of 1/2 for identical particles in the final state. In 
Eq. (I3.24P , A/fc is short for . The a -|- a — )■ a -|- a scattering cross-section due to the 
self-interaction is 

1 11/" d^pz f d^Pi 2 I/O \4e4/ 



^ (3.25) 



647r 



where the last equality holds in the non-relativistic limit. The particle density in physical 
space is 

d^p 



n 



(2vr) 



. (3.26) 



If most states are not occupied, Eqs. (13. 241) . (13. 25 p and (I3.26P imply the usual expression for 
the relaxation rate 

r ~ ^ ~ n (7 5w (3.27) 

where 5v is a measure of the velocity dispersion in the fluid. 

However, in the case of cold dark matter axions, the momentum states that are occupied 

are enormously occupied j2|: Mp ~ 10^^ ( loi^'^Gev ) ^ P ~ Pmax(^) ~ Tl^^^^'i Section V. 
For such a situation, consider the factor 

F{U) = (M + l){Af2 + lW3N'4 - UMU^ + 1) (M + 1) (3.28) 

that appears in Eq. (I3.24p . The first term in F{M) multiplies the contribution to < Mi > 
from a{pi) + a{p2) ^ ^(ps) + 0(^4) transitions whereas the second terms multiplies the 
contribution from a{pi) + 0(^2) ciiPs) + 0(^4). By energy conservation, if the initial 
particle states (e.g. pi and P2) are both highly occupied, the final particle states (ps and p^) 
are almost always highly occupied as well. The RHS of Eq. (I3.24p is dominated by regions 
of momentum space where all four Afc are of order J\f ~ 10^^ The rate for a{pi) + 0(^2) 
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a(p3) + a(P4) transitions, per pair of initial state particles, is multiplied by (A/3 + l)(7V4 + l) ~ 
N"^ compared to the case where most states are empty. However, the rate for a{pi) + a{p2) ^ 
a{p^) + 0(^4) is similarly multiplied by order N"^. The two cancel each other since the terms 
of order M'^ cancel out in F{N'). However, the order terms do not cancel out. The 
upshot is that, although the rates at which quanta enter and leave any one state pi are both 
enhanced factors of TV^, the difference between those two rates, i.e. the relaxation rate, is 
multiplied by one factor of Af: 

V^'^^na5vJ\f . (3.29) 

Using this equation, we show in Section V that cold axions barely thermalize near time 
ti, just after they are produced during the QCD phase transition. However the condition 
T « 5uj that defines the particle kinetic regime is only borderline satisfied at time ti and is 
violated afterward. Shortly after ti the axion fluid enters the condensed regime. Relaxation 
in the condensed regime is discussed in the next subsection. 

We would like to conclude this subsection with remarks about gravitational interactions. 
It has long been re cog nized that gravitational interactions do not fit into the usual discussion 
of thermalization |26|. Indeed the total cross section for gravitational scattering is formally 
infinite, and Eqs. f l3.27p and fl3.29p are clearly invalid. Let's use our derivation of the Boltz- 
mann equation to shed some light on this matter. Eq. (I3.2ip is valid for all interactions 
since it only assumes the ordinary rules of quantum mechanics. However, Eq. f l3.23p is valid 
only in the particle kinetic regime. It does not apply to transitions with very small momen- 
tum, and hence very small energy, transfer since the particle kinetic condition fifj t >> 1 
is violated then. The divergence in the cross-section for gravitational scattering is at small 
momentum transfers (forward scattering), where the assumptions underlying the calculation 
of the cross-section are invalid. Even the notion of scattering is dubious there. We should 
restrict the sum in Eq. fl3.23p to transitions that are in the particle kinetic regime. This 
will cut off the divergence in the total cross section for gravitational scattering but the dif- 
ferential cross section will still be peaked in the forward direction. Forward scattering has 
only limited effectiveness in relaxing a momentum distribution since the particle momenta 
change only a little bit at the time. The cross section appearing in Eqs. f l3.27p and f l3.29p is 
the cross-section for large angle scattering since only large angle scattering has an immediate 
effect on the momentum distribution. The issue does not arise for interactions because 
the differential cross-section is not sharply peaked in the forward direction in that case. We 
conclude that one may use Eq. fl3.27p or (I3.29P with 

^» ~ W 

to estimate the relaxation rate by gravitational scattering in the particle kinetic regime. 



D. The condensed regime 

Consider a huge number of particles occupying a small number K of states. The average 
occupation number of those states that are highly occupied is A/" = ^ >> 1. Assume that 
the energy spread of the highly occupied states is small compared to the evolution rate (to 
be determined) of the system: 5uj « T . For transitions between such closely spaced states, 

e"'^?''* = 1 (3.31) 
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in Eq. f l3.2ip . The first order terms in that equation no longer average to zero. We want to 
estimate the relaxation rate under these conditions. 

It may appear at first that we are embarking on a meaningless quest because the condition 
fl^jt « 1 imphes that the energy difference between highly occupied states is too small for 
anyone to tell whether a particle has made a transition from one highly occupied state to 
another in time t. That is indeed so. However, for the case of interest to us, {Vl^^-t)M >> 1 
and it is therefore meaningful to ask whether of order M particles have made a transition 
between highly occupied states. The latter is the question which is relevant to estimating 
the relaxation rate. 

Let us define ci{t) = ai{t)e^'^'*. The Hamiltonian (13. lip implies 



M 

^lit) = E 2 4c,c,e^^^* . (3.32) 



k,i,j=l 



Let us define further 



ci{t) = Ciit) + diit) (3.33) 

where the di{t), like the ci{t), are annihilation operators satisfying canonical equal time com- 
mutation relations and the Ci{t) are complex c-number functions which satisfy the classical 
equations of motion 

M 

Clit) = E ClC^.e^'' . (3.34) 

For the highly occupied cold axion states, the Ci have magnitude of order y/M. The relax- 
ation rate of the highly condensed cold axions is the inverse of the time scale over which 
those Ci{t) change by order y/JJ . 

Since the sum in Eq. ( I3.34p is dominated by terms for which k, i and j label highly 
occupied axion states, 

^ 1 

k,i,j=l 

For X(j)'^ self-interactions, we substitute Eq. ( 13. 7p . This yields 



P2,P3 

where = Pi +P2 — Ps, and the sum is restricted to the K highly occupied states for which 
P ^ Pmax- We may think of the terms on the RHS of Eq. (I3.36P as steps in a random walk 
in complex space. The magnitude of each step is of order A/" 2 and the number of steps is of 
order K^. Hence 

' Am?V Am?V 

Hence our estimate for the relaxation rate due to Xcf)^ self-interactions in the condensed 
regime [2|: 

Fa ~ ^nXm-^ (3.38) 
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where n = N/V is the density of particles in the highly occupied closely spaced states. 
Likewise, using Eq. f l3.10p . we find that through gravitational self-interactions 



(3.39) 



P2,P3 

The corresponding relaxation rate is 

Tg ~ AnGnm'^f (3.40) 

where i ~ is the correlation length of the particles. 

We note that at the boundary between the particle kinetic and condensed regimes, where 
6uj ~ r, the two estimates of the relaxation rate agree with one another. Indeed, at that 
boundary, 

SvAf^6v-—^ tF ■ (3-41) 

Substituting this into Eq. fl3.29p for the relaxation rate due to Xcp"^ self-interactions in the 
particle kinetic regime yields Eq. fl3.38p which is the corresponding estimate in the condensed 
regime. The same holds true for the gravitational self-interactions. 

Let us also note that Eqs. f l3.38p and fl3.40p are not valid when almost all axions are in a 
single state {K = 1), as when the Bose- Einstein condensation has been completed. Indeed, 
if i^T = 1, there is only one term in the sum on the RHS of Eqs. fl3.36p and fl3.39p that is 
enhanced by large occupation numbers, i.e. the term for which both p2 and ps equal the 
momentum of the single highly occupied state, and it describes an interaction with zero 
momentum transfer. Thus, once the Bose-Einstein condensation is complete and all axions 
are in the lowest energy state, any further thermalization is suppressed. 

Finally consider transitions a(pi) + 0(^2) ^ o^ip^) + aipi) where p2 and are momenta of 
highly occupied states but pi and are not. Such transitions are in the condensed regime 
because the momentum transfer, and hence the energy transfer, is small. Eqs. fl3.36p and 
f l3.39p apply to such transitions and imply that the rate at which states with p > p^aax 
modify their occupation numbers is also given by Eqs. fl3.38p and fl3.40p with the proviso 
that the quanta can only move between states differing in momentum by less than Pmax- 

Eq. (I3.40p has a simple interpretation. The axions, having energy density mn and corre- 
lation length i produce gravitational fields g ~ AirGpi. The rate at which these gravitational 
fields modify the momentum of a particle by an amount Ap is of order 

r„ ~ 4^ ~ AnGnmi^ (3.42) 
^ Ap Ap ^ ^ 

where u is the energy of the particle. For the axions themselves, we obtain the relaxation 
rate Eq. f l3.4Up by substituting uj = m and Ap ~ . Eq. f l3.42p shows that the momentum 
distribution of any particle species is modified by the gravitational fields of the cold axion 
fiuid and therefore that gravitational interactions may produce thermal contact between the 
cold axions and other particle species. This is discussed in the next subsection. 



E. Other species 

Our purpose in this subsection is to estimate the gravitational interaction rates of other 
species - baryons, relativistic axions and photons - with the cold axion fluid. We are moti- 
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vated by the question, discussed in Section V, whether these other species come into thermal 
contact with the cold axion fluid. 

The Hamiltonian describing gravitational interactions between the cold axions and any 
other species has the general form: 



M S 



^=Xl^i4% + J2^'-^rbr + Yl ^^fc«4«!«*«i + ^b'ks^bl^^A , (3.43) 

j=l r=l j,k,r,s 

where A^*^^^ = (A^'*^-^)*. The br are the annihilation operators for quanta of the new species. 
They satisfy canonical commutation or anti-commutation relations. The Ur are the energies 
of those quanta. The other symbols {uj, aj and A^-'J have the same meaning as in Eq. (13. lip . 
As before, we quantize in a box of volume V = with periodic boundary conditions. The 
labels of the new particle states are then r = (n, a), giving their momenta p = and their 
spin a. Their energy is a; = >/p ■ p + where is the mass of the new species. 

We deflne Cj(t) = aj(t)e^'^^* as before, and c^(t) = br{t)e^'^''^. The Heisenberg equations of 
motion for the c^(t) are then: 



-^E^r..4c,4e^^-* (3.44) 



where Q'j'^ = Uk + ujg — ujj — Ur- Because 3-momentum is conserved in each interaction, the 
A^*^^^ have the form: 

K^ks = ~K''ks ^Pk+Ps,Pj+Pr ■ (3.45) 

The important contributions in the sum on the RHS of Eq. fl3.44p are from terms in which 
both j and k label highly occupied cold axion states. Therefore 

K 

c's ~ +^ E KksC:C,c', e^(---)* (3.46) 
k,j=i 

with pr = Ps + Pk — Pj- As before, K is the number of highly occupied cold axion states 
and the Ck are deflned by Eq. (13.330 . Again, the sum on the RHS of Eq. (I3.46P represents 
a random walk in complex space. The number of steps is K"^ and the typical step size is 
A;, Afc', where Af, is the typical value of A^'^^i.^. The rate at which all quanta of the new 
species may move to neighboring states seperated in momentum space by less than 5p ~ | 
is therefore 

Tb,sp ~ K^^\b = \bN (3.47) 

where N = KM is the number of cold axions in volume V . The relaxation rate of the new 
species is then 

r..A.iv|-~A./V^ (3.48) 

where Ap is the momentum dispersion of the new species population. (If the momentum 
dispersion is very different in the initial state than in the flnal state, Ap is the larger of the 
two.) Eq. f l3.48p assumes that the h particles are bosons or non-degenerate fermions. If they 
are degenerate fermions, their relaxation rate is suppressed, relative to Eq. ( I3.48p . by Pauli 
blocking. 
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Also, let us reiterate that when most cold axions are in the lowest energy state, implying 
K = 1, thermalization is suppressed compared to the estimate in Eq. f l3.48p . because there 
is only one term in the sum of Eq. fl3.46p in that case and the momentum transfer vanishes 
for that term. 



1. Baryons 

For non-relativistic species, such as baryons and WIMPs, the term in the Hamiltonian 
that describes gravitational interactions with the cold axions is 

H, = -G [£x£x''-^Mp^ , (3.49) 
J \x — x'\ 

where 

PBix, t) = ^Yl '''^'^ ' (3-50) 

n,n',cr 



and TTj-B is the mass of the non-relativistic particle. This yields 

n-,.(n^.a') " TTZd. ^a' > l<^-OJ-j 



where g = pi — pa is the momentum transfer. Since g ~ £ ^, the B particles have relaxation 
rate 

~ AirGnmmB-: , (3.52) 

APB 

where ApB is their momentum dispersion. Eq. (I3.52p assumes that the B particles are 
bosons or non-degenerate fermions, as is the case for baryons and WIMPs. 



2. Hot axions 



For relativistic species, the term that describes gravitational interactions with the cold 
axion fluid is 

Hr = - [ d^x -KpT''^ (3.53) 



where T°'^{x,t) is the stress-energy-momentum tensor of this species and hap is the pertur- 
bation of the space-time metric caused by the cold axions: 





[x,t) 


^2Gj 








\x — x'\ 


hoki 


[x, t) 


= 






[x, t) 




f a x — ^{^k 









(3.54) 



Note that hf^ = 0. For a scalar field 0(a;) 



Hr = - I d^x hlap a"(/> d^^ . (3.55) 
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After some algebra, Eq. fl3.55p yields 



K = +^7^V^ [^^^^ + - 2 ] ' (3-56) 



where q = pi — ps- The relaxation rate for relativistic scalars through gravitational interac- 
tions with the highly occupied low momentum axion modes is thus of order 



47rGW , (3.57) 



since g ~ £ ^ and Ap ~ u. 



3. Photons 

The term that describes the gravitational interactions of photons with the cold axion 
fluid is 

H-y = -J d'xh^^F'^^F^^ (3.58) 

where is the electromagnetic field strength tensor, and the hap are given by Eqs. (I3.54p 
as before. This yields 

A'tg.k) = + v7^, ■ ■ ^ + ^^"^ ^ '"^^ ■ ^^"^ ^ ■ ' ^^-^^^ 

where e*2 and q are the polarization vectors of the initial and final state photons. We find 
therefore 

~ 47rGnmi (3.60) 

for the relaxation rate of photons. It is the same as for relativistic axions, Eq. ( I3.57p . in 
order of magnitude. 



IV. NUMERICAL SIMULATIONS 

As far as we are aware, there has been no detailed discussion of thermalization in the 
condensed regime prior to this work, presumably because most many body systems ther- 
malize in the opposite particle kinetic regime. The cold axion dark matter fluid may be 
the only physical system for which the condensed regime is the relevant one. Since direct 
observation of cold axions is impossible, we cannot verify the validity of our estimates by 
empirical methods. However, we may use numerical simulation as a check. In this section 
we construct a toy model which thermalizes in the condensed regime but is sufficiently sim- 
ple that it can be simulated numerically. The simulation agrees with our estimate of its 
thermalization rate. 

The Hamiltonian for the systems of interest to us has the general form given in Eq. (13. lip 
which we repeat here for convenience: 



^ 1 
^ = Y1 ^A^i + 5Z 4 ■ (4-1) 
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As before, let us define Cj(t) = aj{t)e^^^^. The Cj{t) satisfy the Heisenberg equations of 
motion 

M 

^lit) = E 2 ' (4-2) 

k,i,j=l 

where = ojk + ^i — oji — Uj. The system is in the condensed regime when a large number 
of quanta are in states whose spread 5u in energy is small compared to the evolution rate of 
the system, so that e^^'^^ = 1 for the dominant terms on the RHS of Eq. (14. 2p . We discussed 
thermalization in the condensed regime in Section III.D. Our estimate for the thermalization 
rate is 

r-v^A^A (4.3) 

where / is the number of interaction terms on the RHS of Eq. (14. 2p between highly occupied 
states, M the average occupation number of those states, and A an average value of the A^^i. 
We may expand any state of the system 

|^(t) >= J2 '^{W,t) e-'^(^-^>)* I {AT} > (4.4) 

m 

in terms of an orthonormal set of Fock space states 

M 

\{Af} >= \Af,,Af,,...,AfM >=l[^={a])^^ |0 > . (4.5) 



j=l V J 

These are the eigenstates of the non-interacting Hamiltonian. The corresponding eigenvalues 
are 

M 

E{{M}) = Y,M,u, . (4.6) 

i=i 

The probability for the system to be in state \{M^ > at time t is |\1'({A/'}, t) 
The time evolution equation 

t^\^{t)>=H\^{t)> (4.7) 



implies 



where 



^iW},t) = \^Ze-^''^^'^{U. + l){U, + l)UuK ^{{ZM],t) (4.8) 



4 

i,j,k,l 



liMp = Afp if pj^i,j,k, I 



= Afp + 1 if p = i or j 

= Mp-1 if p = k or I . (4.9) 

The iZ-^} configuration is the same as the {Af} configuration except that two particles 
are moved, from the states k and / to the states i and j. In obtaining Eqs. (14. 8 p and (14. 9 p 
we assumed that the couplings A^^^ vanish when k = I or i = j. If we do not impose this 
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restriction, Eqs. f l4.8p and (14 .Qp are a little more complicated. However the complications 
do not change the estimate of the thermalization rate, so we ignore them. 

The evolution of the system can be determined by solving Eqs. (14. 8 p numerically. For a 
system in the condensed regime, the important terms on the RHS of that equation all have 

e~*^fci* = 1, so that 

^nm,t) = Yl \^^l^/(M + ^m + m^ln{zm,t) . (4.10) 

i,j,k,l 

Recall that the cold axion fluid is in the condensed regime because the axions occupy rela- 
tively few states, each state is hugely occupied {M ~ 10^^), the particle states have a small 
spread in energy {6ijjT~^ « 1) yet, because 6uj AfT^^ » 1, it is meaningful to ask what 
is the distribution of axions over the closely spaced states. The multiple requirements are 
satisfied only because is very large. But if is large, it is difficult to solve Eqs. (14.101) 
numerically. Indeed, the number of system states (the dimension of Hilbert space) is 
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(N + M-1) , , 

D = - — (4.11) 

iV!(M-l)! ^ ^ 

if is the number of particles and M the number of particle states they occupy. Clearly, 
we will only be able to simulate systems with relatively small and M. (The simulation 
described below has M = 5 and N = 50, hence D = 316 251 and A/" = 10.) Since Af 
cannot be very large in a numerical simulation, how then does one simulate a system in the 
condensed regime? 

We may contrive a toy model in the condensed regime by setting A^^^ = when fi^f^ ^ 0. 
In the example studied below, M = 5 and the states are equally spaced: Uj = jui for j = 
1,2, .. 5. The A^^i all vanish except A2I, A24, A'^l and their complex conjugates Afl, Af^ and 
A2I. The toy Hamiltonian is therefore 

H = Uia\ai + 2ujia2a2 + Sujiala-^ + Auia\a4^ + 5a;ia5a5 + Alla2alaia4^ + All*a\a\a2a3 
+ Allala^aia^ + All*a{ala2a4 + Allala\a2a5 + All*alala3a4 . (4-12) 

The toy model is a far cry from a description of the cold axion fluid but it is in the condensed 
regime in the sense that its evolution is governed by Eq. (I4.10p . and that evolution is non 
trivial. We may ask how quickly does the toy system reach thermal equilibrium starting 
from an arbitrary initial state |A/i, A/2, As, A^, A/5 >. The initial state determines the energy 
E through Eq. (14. 6p . The (microcanonical ensemble) thermal averages of the A/j are 

^,= '4^ . (4.13) 

J2^{E - E{{Af})) 

m 

The N'j were calculated numerically. They differ generally from the Bose-Einstein distribu- 
tion 

^^'^^ = eK-M)/T _ 1 ' (4-14) 

where T is temperature and /i is chemical potential, because the toy model is not in the 
thermodynamic limit. 
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We calculate the quantum mechanical averages 



<A/;(t)>=5^A/;|v&({Ar},t)p (4.15) 

as a function of time and ask how quickly do they reach the thermal averages Mj. 
Fig. 1 shows the result of a particular simulation with N = 50. The initial state is 
|20,5,15,5,5 >, hence E = 120 ui. The thermal averages for that value of E are 
{Af} = (15.25,14.72,9.33,6.16,4.53). The evolution was computed using Eqs. flCTjl and 
(^JS^ with A^^ = All = ^34 = O-l- The time step was At = 10-^ sufficiently small that 
energy is conserved within 1%. Fig. 1 shows the < Mj{t) > as a function of t. Their values 
at t = 2.3 are (15.57, 14.61, 9.42, 6.29, 4.59). There is no doubt that the system thermalizes. 
Eq. f l4.3p estimates the thermalization rate to be F ~ -\/3 x 10 x 0.1 = 1.7, i.e. the thermal- 
ization time 1/F is 0.6 . That estimate is consistent with the time evolution shown in Fig. 
1. We ran simulations for a number of different initial conditions. In each case, the N'j{t) 
approached the appropriate thermal averages on a time scale of order 1/F ~ 0.6. 

It may appear unsurprising (and anti-climactic) that the simulations agree with the ther- 
malization rate that we expected. Let us keep in mind, however, that there have been 
surprises in past simulations of thermalization. Fermi, Pasta and Ulam 29l| performed a 



famous numerical experiment where a set of classical coupled oscillators is evolved in time 
according to a Hamiltonian of the form 

^ ^ 5Z ( 2^ + 2"^^i^? ) + 5Z ^^^^^ ^i^j^kXi . (4.16) 

j=l ^ ^ ijkl 

They found that the system does not thermalize even after a very long time. We repeated 
this relatively simple simulation with M = 10 oscillators and also found that the system 
does not approach thermal equilibrium, even after a very long time. It is reassuring that, 
in contrast, the system of bosonic oscillators behaves in the way we expect. As far as we 
know, ours is the first simulation of the approach to thermal equilibrium of a set of coupled 
quantum oscillators. 



V. AXION COSMOLOGY REVISITED 

In this section, we consider how the late thermalization of axions and other species modi- 
fies the standard cosmological model. It is assumed that the estimates of the thermalization 
rates in Section III.D are correct. We find that cold axions may briefiy thermalize after they 
are first produced during the QCD phase transition as a result of their self-interactions. 
However, this first thermalization era promptly ends and the axions are then decoupled until 
the photon temperature reaches approximately 500 eV. At that time the axions thermalize 
through their gravitational self-interactions. They form a BEC, meaning that most axions 
go to the lowest energy state available to them. The axion correlation length grows to sizes 
of order the horizon. The ratio Ta/H, of the cold axion thermalization rate to the Hubble 
rate, grows then as a{t)~^t^, the axions thermalize on ever shorter time scales relative to 
the age of the universe and the axion state tracks the lowest energy state at all times. We 
find that baryons enter into thermal contact with the cold axions soon after the axion cor- 
relation length has grown to be of order the horizon. It appears possible, although by no 
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means certain, that baryons, photons and axions all reach the same temperature at about 
the time of equality between matter and radiation. In that case and assuming that the 
initial photon chemical potential is negligible, the photon temperature drops by the factor 
(|)4 = 0.9036 compared to the standard cosmological model. We derive the implications 
of this for observations of cosmological parameters, specifically the baryon to photon ratio 
Vbbn at the time of primordial nucleosynthesis and the effective number N^s of neutrino 
thermal degrees of freedom at decoupling. We consider whether neutrinos may also come 
into thermal contact with the cold axions, and what this would imply for //bbn and Nes- 
Finally we show that axion rethermalization by gravitational self-interactions is sufficiently 
fast that axions that are about to fall into a galactic gravitational potential well share their 
angular momenta. By almost all going to the lowest energy state for given total angular 
momentum, the axions acquire net overall rotation, implying V x t? 7^ where v{r) is the 
velocity field of the infalling axions. In contrast, the velocity field of WIMP dark matter 
is irrotational. This provides a means to distinguish axions from WIMPs by observing the 
inner caustics of galactic halos. As mentioned in the Introduction, the evidence for caustic 
rings is consistent with axions and inconsistent with WIMPs. 



A. QCD epoch 

Cold axions are produced when the axion mass turns on during the QCD phase transition. 
The critical time is ti defined by m{ti)ti = 1 where m{t) is the axion mass. At temperatures 
well above 1 GeV, m ~ whereas at temperatures well below 100 MeV, m has its zero 
temperature value. A standard calculation yields [16|] 

where Ti is the photon temperature at ti, and fa is the axion decay constant. How many cold 
axions are produced depends in part on whether inflation occurs before or after the Peccei- 
Quinn phase transition [28]. If before, there are contributions to the cold axion density from 
vacuum realignment, string decay and wall decay. The size of the string decay contribution 
has been the matter of debate and controversy. If after, the contributions from string and 
wall decay are absent and the contribution from vacuum realignment may be accidentally 
suppressed if the initial value of the homogenized axion field happens to lie close to the CP 
conserving minimum. Allowing for these and other uncertainties, we write the cold axion 
number density at time ti as 

nih)=X^ . (5.2) 

X is of order two if there is no inflation after the PQ phase transition and the contribution 
from string decay is of the same order as that from vacuum realignment, X is of order 
ten if there is no inflation after the PQ phase transition and the contribution from string 
decay dominates over that from vacuum realignment (the contribution from wall decay is 
thought to be subdominant always), and X is of order one half the square of the sine of the 
misalignment angle if there is inflation after the PQ phase transition. The misalignment 
angle is a random number, between and 27r, which has the same value in our whole visible 
universe in this case. 
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After ti the number of axions is conserved in the sense that all axion number changing 
processes, such as axion decay to two photons, happen on time scales vastly longer than 
the age of the universe. Eqs. (15. ip and (15. 2p imply then the cold axion number density as a 
function of time 



5 



4 ■10".-/ /. V' (a{h) 



Most of the cold axions are non-relativistic shortly after ti because the axion momenta are 
of order j- at time ti and vary with time as a{t)~^. Hence we may obtain the axion energy 



density today [16|, |28 



7 

a„.||™,«.).0.3x(^)' . (5.4) 
The velocity dispersion of cold axions is 

Ht) ~ 4^ (5.5) 

mil ci[t) 

if each axion remains in whatever state it is in, i.e. if axion interactions are negligible. We 
refer to this case as the limit of decoupled cold axions (see subsection II.C.l). If decoupled, 
the average state occupation number of the cold axions is [2| 

Because A/" is huge, we may expect cold axions to be a Bose-Einstein condensate. Their 
energy dispersion Scu = |m(5f )^ is much smaller than the critical temperature 

m = ^ - 300 GeV Xl (^i^) ^ ^ (5.7) 



C(3) J VlO^^GeVy a{t) 

for BEG. Note that Eq. (15. 7p gives the critical temperature for a relativistic {T^ » m) Bose 
gas. 

Actually, there are four conditions for a gas to form a Bose-Einstein condensate: 1) the 
particles are identical bosons, 2) their number is conserved, 3) their phase space density 
in units of = {2TTh)~^ is of order one or larger (equivalently their average quantum 
state occupation number A/" is of order one or larger), and 4) the particles are in thermal 
equilibrium. The first three conditions are clearly satisfied in the case of cold axions. The 
fourth condition is the critical one since it is not obviously satisfied. Indeed, axions are 
thought to be very weakly interacting. 

Axions are in thermal equilibrium if their relaxation rate Ta is large compared to the 
Hubble expansion rate H. We noted in Section HI that the formula for the relaxation 
rate differs whether the particles in the fluid are in the 'particle kinetic' regime or in the 
'condensed' regime. At time ti, cold axions are at the borderline between the two regimes 
since 

M'^i) = ^Mii)('^^(ii))'~^;;^ = ^(ii) • (5.8) 
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To see whether the axions thermahze by self-interactions, we use Eqs. fl3.25p and f l3.29p 
in the particle kinetic regime, and Eq. f l3.38p in the condensed regime. At ti both estimates 
give the same answer [2|], namely: 

T^ih) ~ H{t^) (5.9) 

indicating that the axions barely thermalize at time ti by self-interactions. After ti 
we must use Eq. fl3.38p since we are in the condensed regime then. It informs us that 
T\{t)/H{t) oc a{t)~^t oc t~2, i.e. that even if axions thermalize at time ti they stop doing so 
shortly thereafter. Nothing much changes as a result of this brief epoch of thermalization 
since in either case, whether it occurs or not, the correlation length i{t) ~ tia{t)/a{ti). 

To see whether gravitational interactions cause the cold axions to form a BEC, we use 
Eq. (ICTj) . It implies 

T„) ^ 8.G,™V* ^ 5 . 10- (^) ' (5.10) 

once the axion mass has reached its zero temperature value, shortly after ti. Gravitational 
self-interactions are too slow to cause thermalization of cold axions near the QCD phase 
transition but, because Tg/H oc a~^(t)t oc a(t), they do cause the cold axions to thermalize 
later on 

B. Axion BEC formation 

The RHS of Eq. (15.101) reaches one at a time teEC when the photon temperature is of 
order 

Te.c~ 500 eVx(^54^)' . (5.11) 

The axions thermalize then and form a BEC as a result of their gravitational self-interactions. 
The whole idea may seem far-fetched because we are used to think that gravitational inter- 
actions among particles are negligible. The axion case is special, however, because almost all 
particles are in a small number of states with very long de Broglie wavelength, and gravity 
is long range. By gravitational self-interactions the axions modify their momentum distri- 
bution till their entropy is maximized for the available energy, which in this case means that 
they form a BEC. 

Axion BEC causes the correlation length to increase. Indeed in an infinite volume, when 
all particles are in the lowest energy state, the momentum dispersion is theoretically zero 
and the correlation length infinite. This ideal state never occurs because thermalization and 
hence BEC formation are constrained by causality. The axions in one horizon are unaware 
of the doings of axions in the next horizon. Hence we expect the correlation length £, which 
may now be thought of as the size of condensate patches, to become of order but less than 
the horizon. The growth in the correlation length causes the thermalization to accelerate; see 
Eq. (15.101) . Once / is some fraction of t, Tg{t)/H{t) cx a(t)~^t^, implying that thermalization 
occurs on ever shorter time scales compared to the Hubble time. 
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C. Thermal contact with other species 



The rate at which non-relativistic species such as baryons change their momentum dis- 
tribution through gravitational interactions with the cold axion fuid is given by Eq. fl3.52p . 
The momentum dispersion of baryons is of order Ap ~ i/Sm^T where T is the photon 
temperature. We will assume here that cold axions are the bulk of the dark matter. The 
Friedmann equation implies then 

1 

for t < teq, where teq is the time of equality between matter and radiation. Hence 



t V'^BEC / V ^1 



'EEC 

Eq. f l5.13p shows that baryons reach thermal contact with the axion BEC when i becomes 
of order t or soon after that. 

Photons are in thermal contact with the baryons, but the nature and degree of this 
thermal contact are changing at the time of axion BEC formation [s^. Baryons interact 
with electrons by Coulomb scattering. Electrons interact with photons by Compton scatter- 
ing, double Compton scattering and bremsstrahlung. Above approximately 1 keV photon 
temperature, double Compton scattering and bremsstrahlung assure chemical equilibrium 
between baryons and photons (the number of photons is not conserved in these processes). 
Below approximately 1 keV photon temperature, Compton scattering is the only important 
interaction remaining. It maintains kinetic, but not chemical, equilibrium between baryons 
and photons till approximately 100 eV photon temperature. Below 100 eV, the degree 
of kinetic equilibrium progressively diminishes till approximately 2 eV, when it disappears 
altogether. 

In any case, as long as there is only thermal contact between baryons and a few low 
momentum modes of the axion field, only a very small amount of energy can be exchanged 
between the axion field and the other species. However, as time goes on, higher and higher 
momentum modes of the axion field reach thermal contact with its highly occupied low 
momentum modes. Eq. fl3.57p gives the relaxation rate of the axion field as a whole, 
including relativistic states. The relaxation rate of photons F^, Eq. (13.601) . is of the same 
order of magnitude. Using Eq. (I5.12p . we have 

r.,/H^TJH^-{^-j (5.14) 

for t < tcq- After equality (t > tcq), ^r/H oc nit oc a{t)^H'^ = constant. Hence, if Tr/H 
does not reach one before equality, it will remain less than one forever afterwards. 

Hence we consider at this stage two possibilities, which we call cases A and B. In case A, 
Tr/H does not reach one before equality, (perhaps because i, although proportional to t is 
much less than t, e.g. i = t/100), and hence thermal contact gets established only between 
baryons and low momentum modes of the axion field. In case B, Tr/H does reach one 
before equality and thermal equilibrium is reached between baryons, axions and photons. 
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This equilibrium is kinetic only since gravitational interactions conserve particle number for 
all the species involved. 

We should ask whether neutrinos may also reach thermal contact with the highly occupied 
low momentum axion modes, in which case neutrinos, axions, baryons and photons would 
all reach the same temperature. We believe this possibility, which we call case C, unlikely 
for the following reason. Eq. f l3.42p does not apply to degenerate fermions because of Pauli 
blocking. Cosmic neutrinos are semi-degenerate since they have a thermal distribution 
with zero chemical potential. Because of partial Pauli blocking, their thermalization is 
slower than that of relativistic axions, Eq. f l3.57p . Since relativistic axions only barely reach 
thermal contact with the highly occupied low momentum modes of the axion field if they 
do so at all, and thermal contact between those low momentum modes of the axion field 
and neutrinos is delayed relative to relativistic axions, it appears most likely that neutrinos 
remain decoupled from the axions at all times. Although we believe case C unlikely, we 
consider its implications below along with the other two cases. 



D. Implications for cosmological parameters 

In case A, the baryons reach the same temperature as the low momentum states of the 
axion field but only a very small amount of energy is transferred from baryons to axions, 
because the relativistic states of the axion field remain decoupled from its low momentum 
states. The cosmological parameters are the same as in the standard cosmological model 
with WIMP dark matter. 

In case B, baryons, axions and photons all reach the same temperature. The axions are 
heated whereas the photons are cooled. The number of photons is unchanged implying that 
their chemical potential stays at zero or, if photons have a small negative chemical potential 
to start with, increases till it reaches zero. There is photon Bose-Einstein condensation, i.e. 
the excess photons condense into the lowest energy available state, which is a plasma wave 
with vanishing wavevector. With zero chemical potential, the photon energy spectrum is 
Planckian, as required by observation. Energy conservation implies 

2 2 

p,. = ^T^^, = p,/ + p,^ = ^4(2 + 1) . (5.15) 
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Hence T^f = = 0.9036 T^, in case B. In Eq. f l5.15p we neglected the contributions 

from the initial state axions, the baryons, and the condensed photons in the final state. 
This is justified because these contributions are respectively of order 10~^^, 10~^ and 10~^ 
relative to the contributions that are included in Eq. fl5.15p . We also assumed that the 
initial chemical potential of the photons is zero. Finally, we ignored the fact that, while the 
transfer of energy between photons and axions takes place, the expansion of the universe 
affects photons and axions differently because axions have mass. The corresponding error is 
of order m/Teq ~ 10~^, or less. 

After the axions have been heated by thermal contact with the photons, the fraction of 
axions in thermally excited states is of order 



V si flO^^ GeV\^ 

4-10-^— . (5.16) 
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All others are condensed in the lowest energy available state. The condensed axions are 
cold dark matter, with the special properties that are the topic of this paper. The axions 
in thermally excited states contribute one bosonic degree of freedom to radiation. After 
axions, baryons and photons have all reached the same temperature, further thermalization 
is suppressed because, once most of the axions are in the same state, the only interactions 
that are enhanced by large occupation numbers are interactions without momentum transfer 
(see Section III.D). 

In case B, cosmological parameters are modified compared to their values in the standard 
model with WIMP dark matter. Since the photons cool between the epoch of primordial 
nucleosynthesis (BBN, for short) and decoupling, and their temperature at decoupling is 

known, their temperature at BBN is larger by the factor (^|) " compared to the standard 
model. Because baryon number is conserved, the baryon to photon ratio rj = n^/n^is smaller 
at BBN than at decoupling. The baryon to photon ratio at decoupling is reliably determined 
from measurements of the anisotropy spectrum of the cosmic microwave radiation. The latest 
WMAP result |3l| is rjdec = (6.190±0.145) 10"^°. Since the photon number density oc 
the value of r] at BBN is 



'7BBN 



Vde 



(4.57 ±0.11) 10 



-10 



(5.17) 



in case B, whereas //bbn = Vdec in the standard model. Under the assumption that there are 
three neutrino species, t/bbn is the main parameter controlling the primordial abundances 
of the light elements. In the standard model, there is generally good agreement between 
the observed and predicted abundances of three light elements (D, ^He, ^He) but there is a 



discrepancy for ^Li, the so-called 'Lithium problem' [32|. The Lithium problem is alleviated 
in case B, and perhaps solved altogether, because of the lower value of r^BBN- On the other 
hand, the agreement that occurs in the standard model between the observed and predicted 
abundance of deuterium is spoiled (lo| . 

Case B differs also from the standard cosmological model in that it has more radiation 
in coUisionless species (axions and neutrinos). The radiation content of the universe is 
commonly given in terms of the effective number Ngg of thermally excited neutrino degrees 
of freedom, defined by 

7/4 



Prad — P-y 
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(5.18) 



where prad is the total energy density in radiation (including photons, neutrinos and axions) 
and p^ is the energy density in photons only. The standard model predicts A^efr = 3.046, 
slightly larger than 3 because the three neutrinos heat up a little during e~^e~ annihilation. 
In case B 



Prad = P7 + Pa + Pi. 



P7 



1 + - + (3.046)- ( — ) - 
2 ^ ^8U1/ 2 
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(5.19) 



taking account of the fact that not only is there an extra species of radiation (axions) but 
also the contribution of neutrinos is boosted because the photons have cooled relative to 
them. Eq. fl5.19p implies N^s = 6.77 . At present, most measurements place A^efr between 4 
and 5, with one a uncertainties of order 1 [31I, The tendency for the measured values 
to be larger than 3.046 has been taken sufficiently seriously to prompt proposals for new 
physics involving extra neutrino species or a neutrino asymmetry 3J]. The Planck mission 
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is expected to measure Ncs with much greater precision 35|. In so doing, it may shed hght 
on the nature of dark matter. 

In case C, the cosmological parameters differ from their standard values even more than 
in case B. Energy conservation imphes in case C 
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T;(2 + l) + p.;(T;,/i,^ 



(5.20) 



where Tj is the common final temperature of photons, axions and neutrinos, and yUj^j is the 
final chemical potential of the neutrinos. In writing Eq. fl5.20p we make the same approx- 
imations as for Eq. (15.151) . When gravitational interactions establish kinetic equilibrium 
between axions, photons and neutrinos, the photons cool whereas the axions and neutrinos 
heat up. Since the number of neutrinos is conserved in this process, the neutrinos acquire a 
negative chemical potential. To obtain Tj we solve numerically 
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with T^i = (^) » T^i, and Eq. fl5:20|l with 

p,f = 2(3.046) 



k 



(27r 



|3 (k-tJ-^f) 
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(5.22) 
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This yields Tj = 0.873T^j and fi,^f = — 0.65Ty. The photons cool even more than in case B 
(Tf = 0.904Tyj) because their heat is transferred to neutrinos as well as axions. In case C, 
the baryon to photon ratio at the epoch or primordial nucleosynthesis is 



Vbbn = (0.873) Vc = (4.12 ± 0.10)10 



-10 



(5.23) 



using again the WMAP value ?7dec = (6.190 ± 0.145)10 The total radiation density is 
also higher. Equating the RHS of Eq. (Km with the RHS of Eq. (Km yields N^s = 8.3. 



E. Tidal torquing with axion BEC 

Let us consider a region of size L inside of which the axion state (i.e. the state that most 
axions are in) stops being the lowest energy available state because the background is time 
dependent. Under what conditions is thermalization by gravitational self- interactions suffi- 
ciently fast that the condensed axions remain in the lowest energy state as the background 
evolves? We use the same heuristic reasoning that led us to Eq. (13.421) . which we later 
verified by more formal derivations in a number of cases in Section III.E. We expect that 
the axion BEC rethermalizes provided the gravitational forces produced by the BEC are 
larger than the typical rate p of change of axion momenta required for the axions to remain 
in the lowest energy state. The gravitational forces are of order AirGnrn^i. In this case, the 
correlation length i must be taken to be of order the size L of the region of interest since 



32 



the gravitational fields due to axion BEC outside the region do not help the thermalization 
of the axions within the region. Hence the condition is 

4nGnm^L>p . (5.24) 

We now apply this criterion to the question whether axions rethermalize sufficiently quickly 
that they share angular momentum when they are about to fall into a galactic gravitational 
potential well. 

We use the self-similar infall model of galactic halo formation to estimate L and p. L is of 
order a few times the turnaround radius R(t), say L(t) ~ 3R(t), whereas p{t) ~ "m-f rot (^) Jmax 
where frot is the rotation velocity and jmax is the dimensionless number characterizing the 
amount of angular momentum of the halo. In the self-similar model, frot(^) ~ R{t)/t and 
R{t) oc ti"*"^ where e is in the range 0.25 to 0.35 0|. Assuming that most of the dark matter 
is axions, the Friedmann equation implies 

AiiGnm = ^H{tf = ^ (5.25) 

for t > teq- The LHS of Eq. fl5.24p is therefore of order 

R(t) , , 1 , , 

2m^ ~ 2mv,ot{t)j (5.26) 

whereas its RHS is of order 

12 1 

= mj^^^v,ot{t)-{— - -) . (5.27) 

The typical value of jmax is 0.18 . Hence Eq. f l5.24p is satisfied at all times from equality till 
today by a margin of order - — A — rr ~ 30 . 

JmaxVg,, 3) 

We conclude that the axion BEC does rethermalize before falling into the gravitational 
potential well of a galaxy. Most axions go to the lowest energy state consistent with the 
total angular momentum acquired from neighboring inhomogeneities through tidal torquing 
[36| . That state is a state of rigid rotation on the turnaround sphere, implying V x w 7^ 
where v is the velocity field of the infalling axions. In contrast, the velocity field of WIMP 
dark matter is irrotational. The inner caustics of galactic halos are different in the two cases. 
Axions produce caustic rings whereas WIMPs produce the 'tent-like' caustics described in 
ref. [5|. There is evidence for the existence of caustic rings in various galaxies at the radii 
predicted by the self-similar infall model. For a review of this evidence see ref. [8|. It is 
shown in ref. ^ that the phase space structure of galactic halos implied by the evidence for 
caustic rings is precisely and in all respects that predicted by the assumption that the dark 
matter is a rethermalizing BEC. 



di 



(to) 



t \ ^ + 9- 



to 



VI. SUMMARY 

The purpose of this paper was to investigate when and to what extent cold axions differ 
from ordinary cold dark matter, such as WIMPs and sterile neutrinos. 

In Section II, we showed that cold axions behave as ordinary cold dark matter on all 
scales of observational interest when they are decoupled. Observable differences between 
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cold axions and ordinary cold dark matter occur only when the axions self-interact. When 
the axions self-interact at a sufficiently high rate, they thermalize and form a Bose-Einstein 
condensate. Before that time, the axions are described by a free classical field and are 
indistinguishable from ordinary cold dark matter on all scales of observational interest. After 
Bose-Einstein condensation, almost all axions are in the same state. If that state is time 
independent, axions again behave as ordinary cold dark matter on all scales of observational 
interest. Observable differences occur if and only if the axions rethermalize so that the axion 
state tracks the lowest energy state. 

In Section III, we calculated the thermahzation rates of cold axions through A0^-type 
self-interactions and through gravitational self-interactions. We described the axion field 
as a set of coupled quantum-mechanical oscillators and asked how quickly such a set of 
oscillators approaches thermal equilibrium. We found that there are two distinct regimes: 
the 'particle kinetic regime' defined by the condition that the energy dispersion is large 
compared to the interaction rate {5oj » F), and the opposite 'condensed regime' where 
Su << r. We derived the Boltzmann equation as an operator statement in the particle 
kinetic regime. However, cold axion thermahzation occurs almost entirely in the condensed 
regime. We derived estimates for the thermahzation rate in the condensed regime, and 
applied the estimates to cold axions and to other species (baryons, photons and relativistic 
axions) that may come into thermal contact with the cold axions. 

In Section IV, we performed numerical simulations to check our estimates of thermahza- 
tion rates in the condensed regime. We constructed a toy model with 50 quanta distributed 
among 5 oscillators. The exact quantum-mechanical evolution of the toy model was solved 
numerically for a variety of initial conditions. It was found that the toy system thermalizes 
and does so on the time scale predicted by the analytical estimate in Section 111. 

In Section V, we investigated what changes the thermahzation of cold axions brings to 
axion cosmology. It was found that cold axions thermalize by gravitational self-interactions 
when the photon temperature is of order 500 eV. When they thermahze, the cold axions form 
a Bose-Einstein condensate. The axion correlation length grows to a size of order but smaller 
than the horizon. Shortly thereafter, baryons come into thermal contact with the axions. As 
time goes on, increasingly higher momentum modes of both the axion field and the photon 
field come into thermal contact with the cold axions. It is possible, but by no means certain, 
that full thermal equilibrium is established among axions, photons and baryons. This would 

imply a drop in the photon temperature by the factor (|) " = 0.9036, at about the time of 
equality between matter and radiation. We derived the implications of this for cosmological 
parameters, specifically the baryon to photon ratio at the time of primordial nucleosynthesis 
and the effective number of neutrinos (a measure of the radiation content of the universe) 
at the time of decoupling. Finally we showed that cold axions thermalize sufficiently fast by 
gravitational self-interactions that the axions about to fall into a galactic gravitational well 
acquire a state of net overall rotation. In contrast, ordinary cold dark matter falls in with 
an irrotational velocity field. The inner caustics of the galactic halo are different in the two 
cases. The occurrence of caustic rings of dark matter in galactic halos is inconsistent with 
ordinary cold dark matter, but consistent with axion BEG. 
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FIG. 1: Time evolution of the quantum mechanical average particle number < Mj >,j = 1..5, 
in each of the five states of the toy model described in the text, starting with the initial state 
(A/i, ...As) = (20, 5, 15, 5, 5). The dots indicate the thermal averages Afj for the corresponding 
total energy. The system reaches thermal equilibrium on a time scale l/F ~ 0.6. 
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